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Rendering participating media requires significant computation, but the ef- 
fect of volumetric scattering is often eventually smooth. This article proposes 
an innovative analysis of absorption and scattering of local light fields in 
the Fourier domain and derives the corresponding set of operators on the 
covariance matrix of the power spectrum of the light field. This analysis 
brings an efficient prediction tool for the behavior of light along a light path 
in participating media. We leverage this analysis to derive proper frequency 
prediction metrics in 3D by combining per-light path information in the 
volume. 

We demonstrate the use of these metrics to significantly improve the 
convergence of a variety of existing methods for the simulation of multiple 
scattering in participating media. First, we propose an efficient computation 
of second derivatives of the fluence, to be used in methods like irradiance 
caching. Second, we derive proper filters and adaptive sample densities 
for image-space adaptive sampling and reconstruction. Third, we propose 
an adaptive sampling for the integration of scattered illumination to the 
camera. Finally, we improve the convergence of progressive photon beams 
by predicting where the radius of light gathering can stop decreasing. Light 
paths in participating media can be very complex. Our key contribution is 
to show that analyzing local light fields in the Fourier domain reveals the 
consistency of illumination in such media and provides a set of simple and 
useful rules to be used to accelerate existing global illumination methods. 


Categories and Subject Descriptors: 1.3.7 [Computer Graphics]: Three- 
Dimensional Graphics and Realism—Raytracing 


General Terms: Theory, Algorithms 
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1. INTRODUCTION 


Rendering participating media is challenging because of the high 
cost of simulating scattering events. But participating media mostly 
blur out details and decrease the contrast of images: some image re- 
gions appear almost locally constant and light beams are practically 
constant in the direction of light propagation. In this article we in- 
troduce a new frequency analysis of local light fields in participating 
media. We show the effect that volumetric scattering has on lowering 
frequency content and contrast. We derive the associated theoret- 
ical framework and provide tools to optimize participating media 
rendering algorithms using the frequency content of light transport. 

Scattering is a long-standing problem in computer graphics where 
arange of techniques, with varying trade-offs between performance 
and accuracy, have been proposed to simulate the interaction of 
light with participating media. Unbiased methods such as path 
tracing [Lafortune and Willems 1993] and Metropolis light trans- 
port [Veach and Guibas 1997] provide accuracy, but often at a 
prohibitive cost. Photon-mapping-based approaches handle partic- 
ipating media [Jensen and Christensen 1998; Knaus and Zwicker 
2011] with different trade-offs. Methods such as Photon Beams 
[Jarosz et al. 201 1a] efficiently simulate low-order scattering, re- 
lying on the accumulation of illumination primitives (e.g., points 
or beams) to compute images. While some approaches exploit the 
lower-frequency nature of lighting in participating media, to our 
knowledge, there is no existing literature on a priori frequency 
analysis of local light fields in volume transport. 

For nonvolumetric surface-based rendering, Durand et al. [2005] 
introduced a frequency analysis of light transport. We extend this 
framework to characterize the behavior, in the Fourier domain, of 
light traveling and scattering inside participating media. Methods 
exist that use the Fourier transform as a global transform operator in 
3D to decouple the frequencies in the scattering equation [Ishimaru 
1997]. Instead, our extension to the frequency analysis framework 
applies to 4D local light fields, along light paths in the medium. 

We build on covariance analysis [Belcour et al. 2013], an 
efficient and practical representation of the covariance matrix of the 
frequency spectrum of the local light field. It was used to accelerate 
the rendering of motion and defocus blur. The covariance matrix 
representation conveys the required information on the Fourier 
transform of the light field at a very small cost, making it tractable 
for path tracing. 
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Fig. 1. 


(b) predicted 3D covariance 


(c) image-space filters 
and camera rays sample density 


(d) reference ( progressive photon beams,4h,44 min) 
inset: equal time comparison (25 min) 


We propose a frequency analysis of light transport in participating media, a broad theoretical tool that allows improvements in a wide variety of 


algorithms. From the predicted 3D covariance of the fluence in the Fourier domain (b), we derive three different sampling metrics in the 3D volume. We present 
multiple example uses of these metrics: to improve image-space adaptive sampling density and reconstruction, we provide sampling density and reconstruction 


filters (c-top); to improve light integration along camera rays, we evaluate a required number of samples along these rays (c-bottom); to improve progressive 
photon beams, we derive the optimal reconstruction radius based on the frequency content ((a) and (d)). In order to ease comparisons, we scaled the covariance 


graphs (b) and increased the luminosity of insets (d). This scene is composed of a pumpkin modeled by user Mundomupa of blendswap.com provided under 
creative common license CC-BY, and of a housefront modeled by Jeremy Birn. 


In this article, we extend the covariance representation to global 
illumination in participating media, including multiple scattering. 
We show that our new formulae for participating media fit nicely in 
the existing framework when the medium is not optically thick (as 
in subsurface scattering). We use the covariance information of the 
local light field spectrum in 4D along light paths to predict the 3D 
covariance of the windowed spectrum of fluence in volumes with 
participating media. We propose four application scenarios where 
this quantity proves useful (see Figure 1). The contributions of this 
article are as follows. 


—wWe give a local analysis of scattering and absorption in the Fourier 
domain along a light path, in heterogeneous participating media. 
The model is compatible with multiple scattering. 


—wWe provide a compact representation of attenuation and scattering 
in the Fourier domain, using covariance matrices. 


—We combine the covariance from many light paths in the medium 
into usable sampling metrics in 3D. 


—Four different computation scenarios are given that benefit from 
our analysis: computation of second derivatives of the fluence, 
image-space adaptive sampling and reconstruction, adaptive sam- 
pling of scattered illumination along rays from the camera, and 
progressive photon beams. 


Note that the term “spectral analysis” usually has multiple mean- 
ings; it is sometimes used to refer to the eigenanalysis of linear 
operators. In this article the term “spectrum” refers to the frequency 
spectrum of the Fourier transform of light fields. 


2. PREVIOUS WORK 


We categorize related work into research on the frequency analy- 
sis of light transport and on the volume rendering of participating 
media. 
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2.1 


In these methods, the Fourier transform is used as a tool for solving 
the scattering equation at once in the entire domain [Duderstadt 
and Martin 1979; Ishimaru 1997]. Some methods use a different 
basis for certain dimensions, such as the Chebychev basis [Kim 
and Moscoso 2003], or spherical harmonics [Dave 1970]. These 
methods in general depend on a combination of very specific con- 
straints: infinite or spherical domains [Dave and Gazdag 1970], 
periodic boundary conditions [Ritchie et al. 1997], isotropic scatter- 
ing functions [Rybicki 1971], and mostly homogeneous scattering 
functions. These conditions make such methods not very suitable 
to computer-generated images where the constraints of uniformity 
and periodicity can hardly be satisfied. 

Our approach is fundamentally different: we use the Fourier trans- 
form as a local tool in the 4D ray space to predict bandwidth—as 
opposed to globally solving the equations—which allows us to han- 
dle nonhomogeneous participating media. 


Fourier Domain Methods for Scattering 


2.2 Volume Rendering 


The field of rendering participating media has a long history. Vol- 
ume rendering based on ray-tracing techniques was first proposed 
for forward path-tracing integration [Kajiya and von Herzen 1984]. 
It has been expanded afterwards to other integration schemes: 
Lafortune and Willems [1996] extended bidirectional path trac- 
ing; Pauly et al. [2000] adapted Metropolis for participating media. 
Photon mapping [Jensen and Christensen 1998] has been shown ef- 
ficient in generating high-frequency light patterns such as caustics. 
Cerezo et al. [2005] survey the state-of-the-art, though the paper is a 
bit dated. 

Recently, various extensions to photon mapping and progres- 
sive photon mapping use photon beams, rather than point sampling 
along rays, to greatly improve the performance of volume render- 
ing [Jarosz et al. 2008b, 201 1a, 2011b; Knaus and Zwicker 2011]. 
These methods, however, remain unaware of image complexity and 
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rely on an accumulation of illumination primitives (e.g., points or 
beams) to compute an image that will eventually be very smooth. 

Several Virtual-Point-Light (VPL)-based algorithms for volumet- 
ric rendering trade off quality and performance. Light cuts and vari- 
ants [Walter et al. 2005, 2006, 2012] achieve scalable rendering 
of complex lighting with many VPLs for motion blur, depth of 
field, and volumetric media (including for oriented media [Jakob 
et al. 2010]). These scalable approaches couple error-bounded ap- 
proximations with simple perceptual metrics. For interactive VPL 
rendering of participating media, Engelhardt et al. [2010] introduce 
a GPU-friendly bias compensation algorithm. Novak et al. [2012a] 
spread the energy of virtual lights along both light and camera rays, 
significantly diminishing noise caused by singularities. 

Multiple approaches aim at efficiently computing low-order scat- 
tering in refractive media. Walter et al. [2009] compute single scat- 
tering in refractive homogeneous media with triangle boundaries. 
Thrke et al. [2007] solve the eikonal equation with wavefront tracing 
for inhomogeneous media with varying refractive indices and Sun 
et al. [2010] develop a line gathering algorithm to integrate com- 
plex multiple reflection/refraction and single-scattering volumetric 
effects for homogeneous media. 


2.3 Efficient Sampling and Reconstruction Methods 


Some works perform adaptive sampling or local filtering using 
heuristics based on the frequency of light transport without explic- 
itly computing frequency information. Adaptive sampling for single 
scattering [Engelhardt and Dachsbacher 2010] permits resampling 
when detecting an occlusion. This approach finds epipolar lines 
and sparsely samples and interpolates along these lines, but finds 
occlusion boundaries to preserve high-frequency details. An epipo- 
lar coordinate system [Baran et al. 2010; Chen et al. 2011] allows 
to interactively compute single scattering in volumetric media by 
exploiting the regularity of the visibility function. 

The structure of the light field [Levoy and Hanrahan 1986; Gortler 
et al. 1986] can be exploited to perform adaptive sampling or 
reconstruction. For surface radiance computation, Lehtinen et al. 
[2011] exploit anisotropy in the temporal light field to efficiently 
reuse samples between pixels, as well as perform visibility-aware 
anisotropic reconstruction to indirect illumination, ambient occlu- 
sion, and glossy reflections. Ramamoorthi et al. [2012] derived a 
theory of Monte Carlo visibility sampling to decide on the best 
sampling strategies depending on a particular geometric configura- 
tion. Mehta et al. [2012] derive the sampling rates and filter sizes 
to reconstruct soft shadows from a theoretical analysis to consider 
axis-aligned filtering. 

Irradiance caching methods [Jarosz et al. 2008a] inherently per- 
form filtering in the space of the irradiance by looking at the ir- 
radiance gradient. For example, Ribardiere et al. [2011] perform 
adaptive irradiance caching for volumetric rendering. They predict 
variations of the irradiance and map an ellipsoid to define the non- 
variation zone with respect to a local frame. 


2.4 Frequency Analysis of Light Transport 


In their frequency analysis of light transport, Durand et al. [2005] 
studied the frequency response of the radiance function to various 
radiative transport phenomena (such as transport, occlusion, and 
reflection). Other works on this subject [Egan et al. 2009; Soler et al. 
2009; Belcour and Soler 2011; Bagher et al. 2012] have enriched 
the number of effects to be studied (motion, lens) and showed that 
filtering and adaptive sampling methods can benefit from frequency 
analysis. Yet, some radiative phenomena have not been studied 
in this framework, including refraction and scattering. We aim to 


Ray in the local 
light field of w 


Central ray 


Fig. 2. Parameterization of a local radiance light field around a ray of 
direction w. We use du, dv as the transverse spatial coordinates of the ray 
and 60, ô¢ as its angular coordinates. 


fill a part of this gap by bringing comprehension of the frequency 
equivalent of volume scattering and attenuation operators, as well 
as showing the usefulness of such analysis with a few practical 
applications. 

A frequency analysis has been carried out for shadows specifi- 
cally by Egan et al. in 4D to build sheared reconstruction filters for 
complex visibility situations [Egan et al. 201 1a], or in the case of 
occlusion by distant illumination [Egan et al. 201 1b]. 


3. BACKGROUND: COVARIANCE OF LOCAL 
SPECTRUM 


Our ultimate goal is to provide a general, efficient tool for predicting 
the variations of the illumination, at different stages of the calcula- 
tion of global illumination, so as to make sampling and reconstruc- 
tion methods the most efficient possible. In prior work [Belcour 
et al. 2013], it was demonstrated that the covariance of the spec- 
trum of the local light field along rays does this job. In this article, 
we perform the mathematical analysis to extend this approach to 
multiple scattering in participating media. This section recalls the 
basics about local light fields, Fourier analysis of light transport, 
and the covariance representation of the spectrum as background. 


3.1 Local Light Fields 


We call the local light field the 4D field of radiance in the 4D 
neighborhood of a ray. Our space of study is the 4D domain of 
tangential positions around a central ray [Igehy 1999; Wand and 
Straßer 2003]. It is parameterized by two spatial and two angular 
coordinates, defined with respect to the plane orthogonal to the ray 
at a 3D position x (see Figure 2). 


3.2 Fourier Analysis 


Durand et al. analyzed the various effects a local light field un- 
dergoes along a light path [Durand et al. 2005]. They showed that 
the effect of light transport operators such as reflection, free space 
transport, and occlusion all correspond to simple operators on the 
Fourier spectrum of the light field. These operators and their equiv- 
alent operator in the Fourier domain are listed in Table I. 


3.3 Covariance Representation 


It is not practical to perform a complete calculation of the full 
spectrum, especially on a per-light path basis. Fortunately, we do 
not need the full spectrum of local light fields to perform useful 
predictions about how the local light field behaves. The relevant 
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Fig. 3. Operations on the light field and their equivalence on the covariance matrix of the light field’s spectra at various stages of light transport (see Belcour 
et al. [2013] for detailed derivations). The middle row shows each operator applied to the covariance © while the bottom row details the constant matrices 
involved. Occlusion adds spatial covariance (O is the 2D spatial covariance of the occluder), BRDF removes angular covariance (Bgg is the angular covariance 
of the BRDF locally to the reflected direction), while rotation and free space travel only cause a reparameterization. 


Table I. Operators on the Light Field Function along a Light Path, 
and their Fourier Equivalent 


Effect Light space operator Fourier operator 

Space travel Angles— space shear Space— angles shear 

Occlusion Product with visibility Convolution by occluder’s 
spectrum 

Projection Scale Inverse scale 

BRDF integral Convolution by BRDF Product with BRDF’s 


spectrum 


information needed is how far and in which directions the spec- 
trum spreads in the Fourier domain. It was demonstrated that the 
covariance matrix of the power spectrum of the local light field is 
sufficient to maintain this information [Belcour et al. 2013]. 

In the most general case, assuming nonstatic scenes, local light 
fields are defined over space, angle, and time, making the light 
field and its power spectrum a 5D function. The effect of motion is 
derived independently of the other 4D operators using a change of 
coordinates in time [Belcour et al. 2013]. In the present document, 
we chose to focus on static scenes only and work with 4D covariance 
matrices. 

For any zero-centered, nonnegative real function f defined over 
the 4D domain, the covariance matrix of f is a 4 x 4 matrix, denoted 
as È and defined by 


1 
Vi j) e{1,..., 4P E= =l (x.e;)(x.e;) f(x)dx. (1) 
Qe xeQ 


In this equation, e; is the i’ h vector of the canonical basis of the 4D 


space Q, while (x.y) is the dot product of vectors x and y. 
The covariance matrix has very interesting properties in terms of 
what we actually need. 


—lIts eigenvectors indicate in which direction function f spreads 
the most and where it spreads the least. 


—lIts eigenvalues are the variance of the function in all four principal 
directions. 


—It is additive, which allows us to accumulate the covariance ma- 
trices of rays. More specifically, the Monte Carlo estimate of a 
composed covariance matrix of many rays is the weighted average 
of the individual covariance matrices with radiance as weights. 


Therefore, the covariance of the power spectrum (amplitude of the 
spectrum) of the local light field can provide us information about 
sampling and integrating the light field function [Belcour et al. 
2013]. In the remaining text, we use the term spectral covariance 
to mean the covariance of the power spectrum of a given quantity. 
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3.4 Relationship with a Gaussian Approximation 


When the covariance © is nondegenerate, it also coincides with the 
covariance matrix of the Gaussian g defined by 
go) = e E, 

Therefore, representing a function by its covariance matrix È is just 
as good as approximating that function by the preceding Gaussian. 
Such an approximation appears very relevant for power spectra, 
which—just like Gaussians—are zero-centered radially symmetric 
functions. 

However, in order to handle degenerate cases, using a covariance 
matrix is more practical than handling Gaussian functions. Besides, 
all transformations we perform over the local light field directly turn 
into algebraic operations on the covariance matrix, with no further 
approximation except for multiplication [Belcour et al. 2013]. 


3.5 Covariance Algebra 


All operators listed in Table I directly act on the covariance matrix 
as algebraic operations, most of which are left-and-right products 
with constant matrices (see Figure 3). Obviously missing in this list 
are the operators to model the effect of volumetric attenuation and 
scattering over the local light field along a light path. We derive them 
in the next section. We also need an efficient way of combining the 
4D light path covariance information in the 3D domain to predict 
how smooth the diffused illumination will eventually be and where 
and how to sample it. This is done in Section 5. 


3.6 Relationship with Second Derivatives 


We show in Appendix A that, for any function f with good regular- 
ity in the neighborhood of a given point xo, the covariance matrix 
of the windowed spectrum of f in the neighborhood of Xo corre- 
sponds to the Hessian matrix of f in the primal domain at Xo, in the 
coordinate system of the principal curvatures: 


1 af (Xo) 
—— | — (x 
4m? f(xo)| | x?” 
This last property will be very useful later on in our analysis as a 
practical method to compute covariance matrices of the windowed 


spectrum of some signals (e.g., the phase functions) around partic- 
ular points. 


Vi Dif) = (2) 


4. FOURIER ANALYSIS FOR PARTICIPATING 
MEDIA 


In this section, we extend the frequency analysis of light transport 
to participating media. We follow a light path, bouncing possibly 
multiple times, into the medium. Our model is therefore compatible 
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Table II. Definitions and Notations Used in the Article 
ôu, ôv, 60, dd Spatial and angular local coordinates 
1(du, dv, 50,5) local light field function 
Qu, Qv, Q9, Qg Spatial and angular variables in Fourier space 
(Qu, Qw, Qa, Qe) Fourier spectrum of the local light field 


x 4D covariance of the local light field spectrum 

Tr 3D covariance of windowed spectrum of fluence 
ôt Coordinate along the central direction of travel 
k(x, y, Z) Volumetric absorption at 3D position (x, y, z) 
Kuv(u, v) Volumetric absorption in plane orthogonal to a ray 
w, Wi, Ws directions of light (general, incident, scattered) 
Pwi, Ws) phase function for directions wi, ws 

p88, 5d) phase function around @;, ws, i.e. P(0, 0) = p(wi, ws) 
Pg Henyey-Greenstein function with parameter g 

a Finite angle for scattered direction 

DUu Convolution operator in the Fourier Q4, &, plane 


Central 
ray 


Fig. 4. Notations for the attenuation operator. We analyze the spectral 
covariance of the attenuation for a small travel distance ôt along the central 
ray. Using small distances allows to assume that the attenuation is constant 
along w. 


with multiple scattering. We derive the equations to model the effect 
of two operators over the frequency spectrum of the local light field 
around this light path: absorption and scattering. We first perform a 
first-order expansion of the phenomena. We then study the Fourier 
equivalent of the two operators and show how to express them using 
algebraic operations on the spectral covariance matrices of the light 
field. Table II summarizes our notations. 

Although the behavior of individual light paths in participating 
media is potentially complicated, we will show that, in the Fourier 
domain, absorption acts like visibility and scattering acts like re- 
flectance. Not only does this fit elegantly into the existing frame- 
work, but it also results in a very simple frequency prediction tool 
for efficiently rendering participating media. 


4.1 Volumetric Absorption 


We first consider the effect of volumetric absorption. When the 
density of particles is not constant in space, energy is not uni- 
formly absorbed as light travels through the medium. This creates 
an increase in spatial frequencies in the signal (similar to shadows), 
which further propagates to the angular domain because of the travel 
of light. We study the effect of volumetric absorption by a density 
function x(x, y, z) acting as an extinction coefficient along a ray, 
for a small travel distance ôt along w (see Figure 4). 

The attenuated light obeys the following differential equa- 
tion [Cerezo et al. 2005]: 


ada 7 = 2) (0) = ria, o). @) 


We perform a first-order approximation of the absorption, con- 
sidering « to be constant for a small distance dt along œ. This allows 
us to integrate this equation as 


U(x + dtw, w) = I(x, w)(1 — 6tk(x)). (4) 


Let Kuv be the restriction of « to the plane, orthogonal to w (which 
means K,,(du, dv) = k(x + duu + dvv)). We adopt the notation 
p(6u, dv) = 1 — dtk,y,(du, dv). In the Fourier domain, Eq. (4) turns 
into a convolution: 


T=T a9, P. (5) 


In this equation, &g,2, denotes a convolution over the spatial com- 
ponent only. The effect of attenuation is therefore identical to oc- 
clusion, except that the mask p = 1 — dfk,, is a function taking 
arbitrary values in [0, 1] instead of a binary function. Let A be the 
covariance matrix of p. Applying the covariance formula for occlu- 
sion (Figure 3), we write the covariance matrix of the convolution 
as the sum of the two covariance matrices. 


D’=D+A (6) 


This equation shows that absorption transfers covariance into the 
spectrum of the local light field. Another way of seeing this is that 
the oscillations of absorption transfer into the light field. 

Computing matrix A in practice is actually simple using Eq. (2): 
we compute the 2D Hessian matrix H(x) of « in the (u, v) basis 
using finite differences and diagonalize it using a 2D rotation R. We 
apply the absolute value and convert it back to the (u, v) coordinate 
system using covariance rotation with the inverse rotation R” (using 
Figure 3): 


5 | [RTIRHCORTRI O 9 
ie 00 0 0 D 
00 00 


It follows that, if a ray crosses a region with transverse sharp tran- 
sitions of the attenuation function (e.g., a transverse transition from 
opaque to nonopaque medium, such as the one depicted on Figure 4), 
the attenuation matrix will represent arbitrarily large frequencies in 
the direction of the discontinuity; this behavior is equivalent to 
binary occlusion. 

Note that, for locally constant and linearly varying volumes, the 
absorption does not affect the spectral covariance of the signal. In 
this case the effect of attenuation is simply the change of the weight 
we will use when combining covariance matrices from multiple 
light paths that we describe in Section 5. 


4.2 Scattering 


In this section we derive the matrix formulation of the change in 
spectral covariance of a local light field along a ray that is scattered 
in a participating medium. Starting from the scattering equation, 
we perform a first-order analysis of the integral and compute the 
Fourier transform of the approximated local light fields. 

The scattering equation used in ray tracing expresses the local 
increase of radiance at x, in direction w, due to light scattering from 
all incoming directions w according to the phase function p [Cerezo 
et al. 2005]: 


IlU + tws, @s)) 
ot 


Integrating for a small traveling distance ôt, we obtain 


(0) = Ks f plo, ws) lx, wdw. (8) 
4r wes? 


otk, 


[x + ôtws, ws) = U(x, Ms) + f p(w, ws) I(x, wdw. (9) 
T Jowes? 
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Fig. 5. Notations for scattering in 3D, with the input and output angular 
coordinate systems aligned. The local light field around the central incoming 
direction w; (respectively scattered direction ws) is parameterized with spa- 
tial and angular coordinates ŝu, ôv, 50, 5@ (respectively du’, dv’, 56’, 5d’). 
The change in spatial coordinates implies a projection of the first spatial 
component only, scaling it by cosa. 


When performing Monte Carlo rendering in participating media, 
the sum on the right is handled by deciding with Russian roulette 
whether the light path scatters or not. Consequently, to study scat- 
tering along a light path that is known to scatter, we need to deal 
with the integral term of the previous equation only. 


4.2.1 Scattering the Local Light Fields. We study the implica- 
tion of this equation in the 4D neighborhood of a couple of directions 
w; and w,, making a finite angle «œ (in other words, cos œ = @;.@;). 
We derive the scattering equation for small perturbations around the 
incoming and outgoing directions. 

We consider the incoming and outgoing light fields to be defined 
in correlated angular frames, for which the first angular component 
lies in the plane containing w, and @;, as depicted in Figure 5. Let 
(60, dd) and (50’, 6d’) be the angular coordinates of the incoming 
and outgoing light fields in these frames and Rsg, sọ the rotation that 
turns the central direction into the local light field direction (50, 5@). 
We also denote by (Su, ôv) (respectively du’, dv’) the spatial com- 
ponents of the 4D frame around w; (respectively ws). 

With this notation in place, we express the scattered contribution 
of the light field around w; to the light field around ws, by restricting 
the integration domain in Eq. (9) to local coordinates around the 
two vectors: 


1,(8u’, Sv’, 86’, 56’) 
otks ; i 
= li(ôu cosa, dv’, 66, dh) p(Rso,3¢@i, Rso',54'Os)- 
36,50 


4r 
(10) 


Note also that in the integrand /; uses outgoing spatial coordinates 
du’, dv’, but incoming angular coordinates 56, 8. Reparameteriz- 
ing spatial coordinates corresponds to a projection of the first spatial 
coordinate, resulting in the 1D scale by cos a that expands the signal 
along the first spatial coordinate. 

We suppose that p(w;, ws) only depends on the relative position 
of the two vectors, in which case we can express it as p(cos œ). 
Given an input angular perturbation (60, 5) and output angular 
perturbation (50’, 5¢’) in the equatorial parametrization, the angle 
a’ between these directions obeys the law of cosines (see, for in- 
stance, Todhunter [1859, page 18] and Figure 5) which, for small 
perturbations, boils down to 


cos(a’) = cos(a + 56; — 80s) cos(5¢; — ês). 
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We adopt the following notation for the phase function in the neigh- 
borhood of (wi, œs): 


p(50’ — 50, 5h’ — 5h) = p(cos(a + 50 — 56’) cos(Sd — 5¢’)). 
Eq. (10) becomes 
1,(6u’, dv’, 0’, 5d’) 


ot 
Zs i 1, (Su! cos a, 5v’, 80, 5b)p(50' — 80, 5’ — 5). 
50,5 


4r 


In the angular domain, this is a 2D convolution between the incident 
light field and the phase function p in the neighborhood of directions 
Wi, Ws. 


4.2.2 Spectrum Covariance Formula. Given this formulation, 
going to Fourier space then follows in a straightforward manner. 
We simply take the Fourier transform on both sides to get 


K, ot r( Qu 
L 


Llu, Qy, Qo, Qg) = 


Qn OE 
4r cosa 


s Qu, Qo, 2) P. 
COS a 


To translate this relationship into covariance matrices, we apply 
the formulae summarized in Figure 3: the convolution adds angu- 
lar bandwidth to the inverse of the covariance (and thus effectively 
removes angular bandwidth), while the scale by 1/ cos œ scales the 
covariance by cos? œ. The outside factor is kept away for normal- 
ization according to Eq. (1) 


Dy = ((Va Ei Va)! + S)`'. (11) 


where § is the covariance matrix of the scattering operator, S de- 
pends on the 2D covariance matrix o, of the windowed Fourier 
spectrum of p, and V, is the scale matrix due to the spatial repa- 
rameterization: 
00 00 cosa 0 
00 00 0 1 
ws 007 ga] Ya=) 0 0 02) 
ool % 0 0 
To compute matrix o,, we use Eq. (2) that directly gives the 
spectral covariance of the phase function around (@;, œs) from the 
Hessian matrix of p (we suppose that the Hessian is diagonal, with- 
out loss of generality; otherwise, an additional rotation is needed 
just like what we did for absorption): 


a2— 


220, 0) 0 


g= S w? 
P ~~ 4?7(0, 0) 0 zoo, 0)| 


(13) 


The effect of scattering is therefore very similar to what a BRDF 
does on the local light field: it removes frequencies. 

It is also interesting to note that, fora = 7, the spectrum co- 
variance in Q, is totally removed by the preceding equation. This 
is because, in this case, the incoming and outgoing directions are 
perpendicular and therefore no variation along u on the incoming 
light affects the outgoing light field. Note also that, for a general 
light path scattering multiple times in a volume, Eq. (11) needs to 
be interleaved with rotations to correctly align coordinate systems 
between two scatter events. 

In summary, we proved that the effect of the scattering operator 
to the covariance matrix will be: a BRDF operator followed by a 
scaling of the spatial component. We will now give an example of 
how to compute S when p is the Henyey-Greenstein function. 
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0 — @ slice of signal 0 — @ slice of signal 0 — @ slice of signal 
Input Scattered, a = 0, g = 0.8 Scattered, œ 0.8 


0 

Predicted covariance 0 5.14 0 0 
of power spectrum 0 0 1.170 0 

0 0 0 1.170 

5.14 0 0 0 

Measured covariance 0 5.14 0 0 
of power spectrum 0 0 0.531 0 
0 0 0 0.531 


Fig. 6. Validation of Eqs. (11) and (14). We measure the covariance of 
the windowed spectrum of an input light beam (left) after scattering with 
two different angles (respectively œ = 0, middle, and a = 0.2, right), 
for a Henyey-Greentein parameter g = 0.8. Top row: Angular slice of the 
signal (in the primal domain). Middle row: Predicted covariance matrix of 
the spectrum. Bottom row: measured covariance matrix, from the 4D data 
sampled with 64+ in [—1, 1]*. Given the order of magnitude of numbers 
involved in the calculation, the measured and predicted values are very 
close (besides, the square root of the diagonals must be compared). Our 
calculation is conservative and well predicts the behavior of the measured 
data. 


4.2.3 Covariance of the Henyey-Greenstein Phase Function. 
There are multiple analytical models of phase functions avail- 
able [Gutierrez et al. 2009]. As a practical example, we give the for- 
mulas of the spectral covariance matrix for the Henyey-Greenstein 
phase function that is most common in the field. This function is de- 
fined using angle a between the incoming and outgoing directions 
as 


1 1- g? 
i 


with cosa = @;.Q,. 


Pea, Os) = 3 


1+ g? —2gcosa)? 


We show in Appendix B that the spectral covariance of the Henyey- 
Greenstein function locally around a, is 


ay l |lhul O 
OME Ge | 0 hzl gi 
with 
h 3g(2Ug? + 1) cosa + g(3 cos(2æ) — 7) 
= 2(g2 — 2g cosa + 1)? 
3g cosa@ 
hy = 


g? —2gcosa+1- 


As expected, the covariance is isotropic for a = 0 (i.e., hy = 
hy) since the Henyey-Greenstein function is rotationally symmetric 
and is null for g = 0, since the function is constant in this case. We 
validate these equations in Figure 6 with a comparison of ground- 
truth and predicted covariance matrices for two different values 
of a. 


Summary. In this section we have derived equations to compute 
the spectral covariance of the local light field along a light path in- 
side participating media. We have shown that absorption increases 
frequency content and acts as the previously defined occlusion oper- 
ator. Scattering low-pass filters the angular frequencies of the input 
local light field with a bandwidth defined by the phase function. 
Thus, it acts like the BRDF operator. 


5. SAMPLING AND INTEGRATION METRICS IN 3D 


In Section 4, we performed an analysis of scattering and attenuation 
in the 4D space of local light fields along rays. In participating 
media, light bounces in all directions and the covariance of a single 
ray cannot be used to predict the overall frequency characteristics 
of the light distribution. In this section we will see how to leverage 
the 4D local analysis to compute a set of sampling metrics in 3D by 
combining the covariance from many rays. 

We consider the following metrics: image-space bandwidth will 
enable efficient image-space sampling and reconstruction; the vari- 
ance in the volume along a ray will prove useful to optimize the 
placement of integration samples for integrating illumination along 
a ray; and finally, a prediction of volumetric bandwidth will predict 
the optimal size of density estimation kernels to improve volumetric 
integration techniques, such as progressive photon mapping, 

Combining the 4D local covariance of many rays into a single 3D 
field also has favorable practical consequences: we favor reusing the 
covariance of a small subset of light paths by storing it in a buffer. 
However, since the spectral covariance of radiance is directional, 
it would ideally need to be stored in a 5D buffer, a potentially 
computationally and memory-intensive datastructure. Fortunately, a 
good compromise is to base our metrics on the frequency covariance 
of the volumetric fluence only, which requires a spatial (3D) buffer, 
at the cost of a very reasonable approximation. 


5.1 The Volumetric Covariance 


The volumetric covariance is the covariance of the power spectrum 
of the fluence in a volume. It bridges the gap between the light path 
formulation of covariance derived in Section 4 and our proposed 
practical improvements of sampling strategies in existing global 
illumination methods. 

We define the volumetric covariance to be the 3 x 3 covariance 
matrix lx, where entry (i, j) is the ij-covariance of the Fourier 
transform of the fluence Fx in the neighborhood of x: 


1 ox 
Doy = —— | Q;2;F(Q)dQ. 15 
Tx) aa | jQ) (15) 
The local fluence F at an offset s around the point x is defined as 


Ris) = | (x +s, d)do. 
wes? 


In practice, the volumetric covariance is computed and stored in 
a voxel grid and the size of the neighborhood considered for each 
voxel is the size of the voxel itself (so, x + s is restricted to lie in 
the voxel). 


Computation. We compute the volumetric covariance by accu- 
mulating contributions of individual light paths traversing the vol- 
ume. At a point x, the 4D spectral covariance È of an incident light 
path in direction w carries the illumination from a very localized set 
of directions around w. The 2D spatial submatrix of the 4D covari- 
ance of the local light field around is therefore a slice of the 3D 
covariance of the integrated radiance in the plane orthogonal to w. 

Consequently, we compute the covariance of the fluence at x by 
summing up the 2D spatial slices of the covariance matrices of each 
incident light path, padded to 3D with zeroes, and rotated to match 
the world coordinate system. Since X lives in Fourier space and 
the summation happens in the primal domain, submatrices need 
to be extracted from =~! and inverted back to Fourier space after 
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Shear + Attenuation 


Shear + Attenuation 


Fig. 7. To estimate the 2D spatial covariance matrix in image space, we 
accumulate slices of the volumetric covariance T along the ray using Eq. (17) 
and the equations of attenuation (6) and (7). 


summation: 


-1 
T,= (/ RE na RUM) . (16) 
wes? 


In this equation, the notation X~!|sx sy refers to the 2D spatial 
submatrix of matrix =~!, while R, is the 3 x 2 matrix converting 
the two local spatial coordinates around w into the three coordinates 
of the world. Finally, 7 (w) is the normalized incident energy along 
incoming direction w. 

In practice, the integral in Eq. (16) is computed using a classical 
Monte Carlo summation, as light paths in the volume cross voxels 
to which they contribute. We do not need to explicitly compute 
I(@) since it is naturally handled by the photon-tracing approach: 
the number of paths crossing a voxel is proportional to the fluence. 
We only record how many times each voxel was hit, for proper 
normalization. 


5.2 |Image-Space Covariance 


We want to compute image-space covariances for adaptive sam- 
pling. The angular submatrix of the covariance © at the camera can 
be used to derive sampling densities and reconstruction filters for 
ray tracing at each pixel [Belcour et al. 2013]. 

The most straightforward method to obtain © for each pixel in 
the screen would be to accumulate covariance matrices from light 
paths reaching the camera, applying the theory of Section 4. While 
this eventually provides an unbiased estimate of the image-space 
covariance, it needs many light paths to obtain a reasonably noise- 
free estimate. 

It is the spatial variation of “light intensity” in the participating 
medium that will show up as angular bandwidth at the camera and, 
after projection, as spatial bandwidth on the screen [Durand et al. 
2005]. Consequently, we propose computing screen-space band- 
width from I’. To compute È at the camera, we slice the volumetric 
covariance I’ orthogonally to camera rays, pad it to 4D with null an- 
gular covariance, and accumulate it using Eqs. (6) and (7) to account 
for the attenuation between points along the ray and the camera. At 
pixel p, corresponding to a ray with origin c and direction d, we 
have the following. 


XU(p) = KOTI (AQ) + Vestalxy)Tadt (17) 


1 D 
JS K@dt Í 


The normalization constant K (t) accounts for how much energy 
is associated with each voxel in the covariance grid. The result- 
ing matrix is a 4D covariance from which we extract the angular 
component at the camera. The process is illustrated in Figure 7. 
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Eq. (17) is an approximation because it implicitly supposes that 
the last bounce of light before the camera has no angular covariance, 
meaning that the last scattering step is isotropic. In practice we found 
this approximation to have no visible effect. 


5.3 Ray-Space Variance 


When gathering energy in participating media, one needs to inte- 
grate illumination along rays. For each ray, we need to make sure 
that the spacing between integration samples avoids aliasing with 
respect to how much the illumination varies in the volume. This 
requires the variance v(t) of the fluence function along the ray. The 
variance of the fluence in a particular direction œw is naturally given 
by a 1D slice of the volumetric covariance matrix, in direction w 
(supposedly a unit vector). 


v(x, w) = wo Tyo (18) 


5.4 Optimal Kernel Size for Density Estimation 


Methods based on density estimation such as photon mapping need 
to carefully adapt the size of the kernel to collect photons: too small 
a kernel increases variance, while too large a kernel increases bias. 
Silverman gives an estimate of the mean integrated square error for 
density estimation (see Silverman’s monograph [1986, page 85]) 
that depends on the Laplacian of the estimated function. In the 
case of photon mapping methods, the function is the fluence and 
Silverman’s estimate gives 


1 
bias, (x) = sah F(x) 


In this formula, h is the radius of the kernel used for density estima- 
tion and a is a constant that only depends on the shape of the kernel. 
We use again Eq. (2) that links the diagonal of the covariance matrix 
to the absolute value of the partial second derivatives of F(x), to 
obtain an upper bound on the absolute value of the Laplacian from 
the trace of the covariance matrix. 


|AF(x)| < 427 FWT) 


Equality holds when the second derivatives share the same sign. 
From this, we have an upper bound on the density estimation bias. 


[bias,,(x)| < 277h? a F(x) tr(Tx) (19) 


This equation directly gives us the largest possible kernel radius h 
to always keep the bias below a given threshold. 


Summary. In this section we have explained how to combine 
the 4D spectral covariance of many rays into volumetric 3D covari- 
ance. We have derived interesting sampling metrics from volumetric 
covariance. In the next section, we will explain how to use these 
metrics to improve existing rendering algorithms in three different 
application scenarios. 


6. IMPROVEMENT OF EXISTING SAMPLING 
AND RECONSTRUCTION METHODS 


In this section we demonstrate the usefulness of our analysis from 
Section 4, as well as the sampling prediction metrics we derived 
in Section 5. We examine four different calculation steps that are 
involved in computational methods of global illumination in partic- 
ipating media and we show that our sampling metrics can be used 
to accelerate them. 
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Fig.8. The volumetric covariance T is the covariance of the fluence around 
each point in the volume. We show values of T from the covariance grid as 
an ellipsoid (iso-value of the 3D Gaussian with same covariance). Top: in 
the region of diffusion, I is really small. Middle: on the shadow boundary, 
T is large along the normal to the shadow volume. Bottom: in the caustic, IT 
is large orthogonally to the direction of the caustic (caution: all graphs are 
normalized to help demonstrate the shape of the covariance. To compare the 
covariance quantitatively, look at the eigenvalues listed for each dimension). 


Our sampling metrics need the computation of volumetric 3D co- 
variance, as defined in Section 5. To compute and store volumetric 
covariance, we use a voxel grid, the covariance grid. In the use cases 
that follows, we always read the values of T in that grid to compute 
the required metrics. All results are computed on an Intel 17-3820 
with four cores at 3.60 GHz per core and an NVidia GeForce GTX 
680. We use eight threads to benefit from hyperthreading. Unless 
noted, we use a 64° covariance grid. The covariance grid population 
algorithms run ona single thread, while we use multiprocessor capa- 
bilities for volume integration (Section 6.3 and 6.4 using OpenMP) 
and for density estimation (Section 6.5 using CUDA). 


6.1 The Covariance Grid 


We sample light paths and populate the covariance grid using 
Eq. (16). We also record how many times each voxel in the grid 
is visited by light paths, so as to maintain information for proper 
normalization. 

This calculation is not view dependent. Depending on the appli- 
cation, we populate the covariance grid using a fixed proportion of 
the light paths used for the simulation (in Section 6.5), or fill it up 
once before the simulation (Sections 6.3 and 6.4). Figure 9 refer- 
ences the values used for the different scenes. For the algorithms of 
Sections 6.3 and 6.4, ray marching and filling the covariance grid 
with 100,000 light paths takes 21 seconds for a 64° covariance grid 
with the Halloween scene. We used as many as 10,000 light paths 
for the Sibenik scene, as the lights are spot lights. With this amount 
of light paths, it took 8 seconds for ray marching and filling for the 
64? covariance grid. 

Figure 8 shows the volumetric covariance predicted by our system 
in three different locations of a scene showing volumetric caustics 
and shadows. 


Halloween Sibenik 
Imane Space Cov. grid samples 100 000 10 000 
ee SP Step size 0.01 0.002 
Cov. grid samples 100 000 10 000 
Pye ram Step sizes 0.1- 0.01 | 0.1 — 0.002 
Naive Step size 0.01 0.002 
Qa 0.7 0.7 
Prog. Ph Beams 
Ha none Samples per pass 5 000 5 000 


Fig. 9. We list the different parameters used for our results section. We 
report the number of light paths used to fill the covariance grid, the distance 
between samples along the eye ray for the integration in the volume, and 
Progressive Photon Beams [Jarosz et al. 201 1b] parameters (radius reduction 
ratio œ, and the number of light paths sampled per pass). We report only 
scenes that are common to the three algorithms. We used os = 0.1, oa = 0.1, 
and g = 0.5 for the parameters of the volume. For our adaptive sampling and 
image-space filtering algorithms, we used 8 jittered supersampling samples 
per pixel to obtain anti-aliasing. 


T T T T 
Fluence 2nd derivative in primal space (3-points estimate) 
Fluence 2nd derivative in primal space (21^3-sized filter) 
Fluence 2nd derivative using 3D covariance (Equation 2) —— | 


E of 


Total rays cast 
2 f f i 1 
o 10000 20000 30000 40000 50000 


Fig. 10. Comparative estimate of the second derivative of the fluence in the 
primal domain (red and green curves) versus Eq. (2) (blue curve), for point 
x in the direction of the arrow, as a function of the number of beams cast. 
In both cases we used grids of 128? voxels to store the covariance and the 
fluence. In the primal space, the second derivative is estimated in two ways 
from the fluence grid. Red curve: using a finite difference scheme between 
immediate neighbor voxels. Green curve: using a very large Gaussian filter 
around the measured voxel (averaging the nearby 21° voxels). For the blue 
curve, we simply applied Eq. (2) to the covariance matrix picked at that 
voxel in the covariance grid, without filtering. With as low as 5,000 total 
rays cast in the scene, our estimate outperforms the costly filtered estimate 
in the primal domain. 


6.2 Efficient Prediction of the Hessian of the Fluence 


The covariance grid naturally enables a stable computation of sec- 
ond derivatives of the fluence. In this section we study the benefit 
of estimating second derivatives from the covariance in frequency 
space using Eq. (2), rather than trying to estimate second derivatives 
in the primal domain using a finite differences scheme. 

We present such a comparison in Figure 10 in a simple volume 
containing a point light source, as well as for two different integra- 
tion schemes: a 3-points second derivative estimate that naturally 
proves very noisy, as well as a heavily filtered estimate using the 
second derivative of a Gaussian over 21° neighbor voxels. This ex- 
periment not only proves that our model gives a correct estimate of 
the second derivatives, but also that it converges faster than meth- 
ods based on the primal domain. This is not surprising, because our 
method does not require explicit differentiation over the path-traced 
illumination the way the primal domain estimate does. 

Methods that rely on linear interpolation between 3D illumination 
samples in the volume theoretically have an interpolation error that 
is proportional to the second derivatives of the illumination. This has 
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been proven for surface irradiance caching methods [Schwarzhaupt 
et al. 2012]. Such an error estimate outperforms existing heuris- 
tics [Jarosz et al. 2008a]. 

Although extending our work to volumetric irradiance caching 
is beyond the scope of this article, we believe that the irradiance 
caching error estimate based on the Hessian can be extended to 
volumes with an error estimate proportional to the Laplacian of the 
fluence (see, for instance, Schwarzhaupt et al. [2012, Eqs. (5) and 
(6)]). Moreover, the shape of the influence regions of irradiance 
cache samples will be ellipsoids aligned with the eigenvectors of 
the Hessian of the fluence in the volume. This opens interesting 
research avenues toward the replacement of existing heuristics for 
error and influence regions for volumes [Ribardiére et al. 2011], 
especially removing the need for special treatment of occlusion 
[Schwarzhaupt et al. 2012] that the covariance analysis naturally 
handles. 


6.3 Image-Adaptive Sampling and Reconstruction 


An effective method for rendering images with varying local band- 
width is to compute image-space converged illumination samples 
and filter these samples using an appropriate 2D reconstruction 
kernel. For an optimal result, it is necessary to know in advance 
the optimal sampling density (in samples per square pixels) and 
the shape of the reconstruction filter at each pixel [Belcour et al. 
2013], or to compute it from a subset of the light paths used in the 
image [Overbeck et al. 2009; Rousselle et al. 2011]. 

The optimal sampling densities N(p) and the shape of the 2D 
image reconstruction filter f, can be derived for each pixel p from 
the covariance &(p) of the local light field at that particular pixel 
[Belcour et al. 2013]: 


N(p)=kyi=(p)| and f(x) = eA EOD, (20) 


In this expression, X(p)~'|x, is the spatial slice of the inverse of 
the spectrum covariance matrix in the image plane at pixel p. In 
other words, &(p) is the covariance of the Gaussian whose variance 
matches the bandwidth of the image according to the Nyquist rate. In 
practice, for each pixel, we trace a single ray through the covariance 
grid and apply Eq. (17) to compute X(p). 

The number of samples per square pixel is proportional to the 
determinant of the screen-space spectrum covariance. Furthermore, 
the shape of the filter is obtained by slicing the covariance of the 
signal along the spatial dimensions at each pixel. The constant k 
lets us express N(p) as a fraction of the total number of samples 
allocated for the entire image. 

To compute the image, we obtain the required number of sam- 
ples per pixel using Eq. (20) and form an image sampling density 
map. We use this map as a probability density to draw pixels to 
be computed, using rejection sampling. For each pixel to compute, 
we estimate the radiance using path tracing. Each image sample is 
therefore a converged illumination value. Finally, we reconstruct the 
image by filtering the image samples around each pixel p with the 
filter f,, that is given by Eq. (20). 

This computation is efficient because it samples the image very 
sparsely when the resulting image is predicted to be smooth. 
Figure 11 shows the result of such a computation for the pumpkin 
scene. The number of computed pixels is 43% of the total number 
of pixels in the image. The results also show that we correctly pre- 
dict the shape and size the of reconstruction filters. For more trivial 
scenes, the gain is even better. 
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Fig. 11. We demonstrate the use of our prediction metrics for image-space 
filtering and reconstruction of single scattering. We predict the varying 
density of image samples to compute using Eq. (20) (/eft) as well as the 
shape of the Gaussian filter to reconstruct from these samples at each pixel 
(middle), to reconstruct the image (right). 
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Fig. 12. Computing the amount of light reaching the camera through a 
given pixel requires integrating the radiance scattered by the participating 
medium towards the camera. Using an arbitrary number of uniformly spaced 
integration samples (red) might do a bad job if not accounting for the 
distribution of light into the volume. We propose instead to distribute samples 
according to the local variance of the fluence along the ray (green) in order 
to adaptively sample high-frequency regions. 


6.4 Adaptive Sampling along a Ray 


For each pixel that we compute, we need to integrate the illumina- 
tion that is scattered towards the camera. For this we integrate the 
radiance after a last scattering event. Instead of using uniform sam- 
ples, we adapt the sampling density to the variance of the fluence 
along the ray to the camera, as computed in Section 5.3. The vari- 
ance is directly computed from the covariance grid using Eq. (18). 
This allows us to place samples in regions where the path crosses 
rapid changes in illumination, which can be caused by volumetric 
shadows, for example, as illustrated in Figure 12. In practice we first 
uniformly sample the radiance along the eye path and then we use 
additional samples in proportion to the local variance estimate along 
the path. 

We provide a simplified algorithm for the adaptive integration 
used with Algorithm 1. It uses Eq. (18) to evaluate the required 
distance between two samples in the volume. 

Since our algorithm takes advantage of adaptive sampling on 
shadow boundaries, we are able to reduce the aliasing of light shafts 
caused by undersampling high-frequency regions. Figure 13 shows 
that our results on the Sibenik cathedral model outperform uniform 
sampling, at equal computation time, in the detailed shafts. 

We summarize the timings of the image-space-adaptive sam- 
pling and the eye-path-adaptive sampling algorithm compared with 
an equal-quality naive ray-tracing approach in Figure 14. Both al- 
gorithms save computation time by adapting the workload to high- 
frequency regions. 
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Fig. 13. We compare our variance-driven integration method to naive uni- 
form sampling, at equal computation time (14 min). Our adaptive sampling 
clearly removes aliasing caused by the shaft from the Rose window. Inset: 
total number of samples used per pixel for our algorithm. (Model of the 
Sibenik cathedral by Marko Dabrovic.) 


scene Image space | Eye space Naive 
Sibenik 19m 14m 1h 40m 
Halloween 7m 6m 30s 22m 


Fig. 14. Our adaptive sampling and image-space filtering approaches save 
computation time compared to a naive ray-tracing approach for the same 
quality. Eye-path-adaptive sampling and the naive implementation use 8 
samples per pixel for anti-aliasing. The image-space method adapts the 
number of samples up to this limit. 


ALGORITHM 1: Our adaptive sampling algorithm computes the 
single-scattering radiance for an eye ray defined by its position x in 
space and direction w. It returns the light scattered by the volume 
in the interval [0, Tmax] along the ray. Our variance estimate v 
(Eq. (18)) provides a step size in the volume. Note that the last 
step requires special treatment as the step might be larger than the 
remaining integration distance. We do not show it in this example 
to keep a compact formulation. 


function ADAPTIVEINTEGRATION(X, w, Tnax, SteP mins SteP max) 
rad=0 
t=0 
while ¢ € [0..T,,,,] do 
Xx =x+ É 


ar 2./v(X;, @) 
dt = clamp(dt, SteP min» SteP max) 
rad = rad + dt integrateRadiance(x,, —w) 
t=t+dt 
end while 
return rad 
end function 


We investigated different resolutions for the covariance grid 
(Figure 15). A size of 64° for the grid was sufficient for all our 
scenes. Coarser grids will provide a conservative estimation of 
frequencies and lead to poor performances, while finer grids will 
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Fig. 15. We analyse the impact of various resolutions of the covariance grid 
(323, 643, and 128°) on the prediction of the required number of samples 
along camera rays. Smaller grid sizes bring more conservative results while 
larger grids are more costly to handle. A size of 64° performs well for most 
scenes. 


Photon 
beam 


Fig. 16. Given a camera ray (green) and a beam, we use the radius rg , esti- 
mated by the covariance analysis, instead of the radius r; of the progressive 
photon mapping, when r; is smaller. Using this, we gather more beams in 
low-frequency regions and decrease the variance of the image. 


naturally increase the cost of ray marching in this structure. The 
cost of filling the grid is linear with the grid-edge size. For the 
Sibenik scene using 10K light paths, ray marching and filling took 
4s for a 323 grid, 8s for a 64° grid, and 17s for a 128° grid. We 
found that a 64? grid provides a good trade-off between construction 
time, quality, and time required to ray-march during rendering (see 
Figure 15) in all our tests, except for San Miguel where we needed a 
256° grid. 


6.5 Improved Multiple-Scattering Photon Beams 


We study the possible improvement of the convergence rate of pro- 
gressive photon beams (PPB) to illustrate the benefits of frequency 
analysis. 

In the original algorithm, photons are traced in the scene con- 
taining a participating medium and the paths of propagation (called 
beams) are stored [Jarosz et al. 2011b]. Then, for each pixel, rays are 
shot from the camera and the density of beams along the ray is esti- 
mated using a 2D circular kernel. This is repeated while decreasing 
kernel size until convergence is satisfactory. 

Just like any other density estimation technique, the PPB algo- 
rithm fights between too small a reconstruction kernel—causing 
variance—and too large a reconstruction kernel—causing bias. 
Whereas the original algorithm keeps reducing the kernel size 
as more beams come along, we can afford to stop reducing it as 
soon as this size ensures that the density estimation bias is be- 
low a certain threshold. We know exactly when this happens from 
Eq. (19). 

During the gathering pass, for each eye ray, we test for its distance 
d to the beams stored (Figure 16). At the closest point to each beam 
along the ray, we look into the covariance matrix and estimate the 
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(b) our improved PPB, 9.8M beams, 3h 10min 


(a) standard PPB, 10M beams, 3h10min 


(c) reference (133M beams) 


Fig. 17. The full San Miguel scene with a rather indirect illumination shows a very smooth volumetric environment in a significantly complex geometry 
(2.5M polygons). In this context, our frequency prediction method makes the progressive photon beams converge much faster. Because of the size of this scene, 
we used 256°-wide covariance and occlusion grids. The most severe limitation of covariance tracing in this scene was due to the occlusion grid size. Too small 
a grid would cause overestimation of visibility and conservatively populate the covariance grid with high frequencies. We used Morgan McGuire export of the 
San Miguel scene. As he noted on his Web site (http: //graphics.cs.williams.edu/data), some geometry (chairs) and textures (walls and columns) 


are missing. 


ideal gathering radius r, using Eq. (19). We gather that beam only 
if 
d < max(/j, ro), 


where r; is the radius given by the photon beam method for pass #i. 
In other words, we replace the gathering radius of progressive pho- 
ton mapping by a specific radius for each pair (eye ray, photon 
beam) that is adapted to the local variations of the signal. This 
adaptive radius computation stops us from decreasing the radius in 
regions of low bandwidth and therefore significantly reduces vari- 
ance while keeping the error uniformly controlled. We implemented 
this gathering in a CUDA kernel. 

We follow a pure ray-tracing approach in our PPB implementa- 
tion. We generate light paths and eye paths on the CPU and transfer 
them on to the GPU. This lets us simplify the generation of rays, to 
follow specular paths from the camera, and to avoid duplicating the 
scene data on the GPU. This explains the timing differences between 
our PPB implementation and the one of Jarosz et al. [201 1b]. 

We validate our improvement of progressive photon beams using 
the San Miguel scene (Figure 17), Halloween scene (Figure 18), 
Cornell box (Figure 19), and Sibenik cathedral (Figure 20) scenes. 
In all cases, our covariance framework correctly estimates the high- 
frequency regions due to the illumination. San Miguel, Halloween, 
and Sibenik (Figures 17, 18, and 20) scenes have significant indi- 
rect illumination; they prove that our method converges faster than 
PPB. San Miguel also demonstrates the scalability of our technique. 
The nonhomogeneous Cornell box scene (Figure 19) validates that 
our analysis and filtering methods correctly handle nonconstant 
scattering parameters. In this example, the scattering parameters 
are varied based on Perlin noise. Image resolutions are reported in 
Table II. 

At equal computation time, we achieve a much better convergence 
in smoother regions of the image while we approximately keep the 
same convergence in the highest-frequency regions such as shaft 
borders and caustics, as predicted. However, our method will always 
eventually stop reducing the kernel size, contrary to the classical 
photon beam approach. It just happens later for higher-frequency 
regions of the image. 
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(a) our improved PPB, 9.7M beams; (b) standard PPB, 10M beams, 25 min 
25 min 


Fig. 18. The Halloween scene combines high-frequency white shafts with a 
low-frequency orange multiple scattering. Our covariance prediction allows 
to filter out the noise due to the diffuse component while preserving edges 
of shafts. 


Improved PPB 


Fig. 19. In this figure, a nonhomogeneous volume is illuminated. The walls 
of the Cornell box provide indirect color bleeding from diffuse reflections. 
We compare our algorithm after a run of 31 min (8.5M beams) with an 
equal-time run of Progressive Photon Beams (10M beams). Our algorithm 
reduces noise thanks to our adaptive density estimation. 
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our method (20M beams, 2h35 min) 


equal-time PPB (20.5M beams) 
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our method 


equal-time PPB 


equal-time PPB 


p JA 


our method 


Fig. 20. Using our analysis metric our algorithm predicts the low-frequency part of the volume where the progressive photon beam can afford to keep large 
collection radii while controlling the bias. As a result, we provide a much better converged image at equal computation time as classical progressive photon 


beams; here, in a scene with multiple scattering. 


Table II. Image Resolutions Used for Our Rendering Tests 


scene San Miguel | Cornell Box Sibenik Halloween 


resolution | 1000 x 1024 512 x 512 714 x 968 | 512 x 512 


7. IMPLEMENTATION AND DISCUSSION 
7.1 Implementation 


Taking advantage of symmetry, 4D covariance matrices only need 
10 floats to store (instead of 16) into light rays equipped with 
covariance information. An additional float is also needed to keep 
track of light intensity for proper normalization, since our definition 
of covariance has no units (see Eq. (1)). The combined covariance 
matrix £2 of two different paths with covariances X; and X> and 
intensities J; and J, after proper alignment, is 


X12 (LXi + hX). 


h+h 
We also store, along with the covariance information, the tangent 
frame of the local parametrization (two 3D normalized vectors). 

A photon that carries covariance is initialized at the light source 
and covariance is updated as the photon path is sampled [Belcour 
et al. 2013]. For instance, a square diffuse light source will pro- 
duce rays with null angular covariance and spatial covariance that 
depends on the size of the square. 

During light propagation, the covariance matrix © of rays is 
updated when the light is reflected, refracted, or occluded using 
the equations of Figure 3 (see Algorithm 2). Each time a ray is 
scattered, we transform its covariance using Eq. (11) and, for each 
absorption event, we apply Eqs. (6) and (7). Eventually, the various 
operators involved boil down to sums, products, and inverses (or 
pseudo-inverses) of 4D covariance matrices, which is carried out 
very efficiently using explicit formulas. 

The covariance grid uses 6-float arrays to store the individual 3D 
spatial covariances I’, and an additional float for the light intensity 
normalization factor per voxel. We populate the covariance grid us- 
ing Eq. (16), basically summing up matrices multiplied by intensity. 


ALGORITHM 2: The tracing of frequency photons is straightfor- 
ward to implement. It only requires that we update the covariance 
matrix at specific steps. Note that it requires the ray-tracing engine to 
compute specific information for intersection with geometry (such 
as local curvature). The R matrix is the factorized matrix of projec- 
tion, alignment, and curvature, both before and after reflection. T 
is the covariance of the texture matrix. 
function TRACEFREQUENCYPHOTON 
{p, œ} < sampleLight() 
x < computeLightCovariance() 
while russianRoulette() do 
p < traceRay(p, œw) 
for all voxels v until hit do 
updateVoxelCovariance(v, X) 


E < TIET 
x<D+0O 
x<D+A 
end for 
@ < sampleBRDF() 
E < RIER; 
EE+T 
SRS + B)`'R, 
end while 


end function 


Depending on the application, the covariance grid is either populated 
once, using a fixed number of rays (for, e.g., adaptive reconstruction 
and ray-space integration) or updated during the computation (for 
progressive photon beams, where only 10% of the light paths carry 
covariance information). 

We did not use ray differentials in our implementation of PPB 
(nor in our improvement). Jarosz et al. [2011b] showed that using 
ray differentials was beneficial for the convergence of specular paths 
from the light. But they use a constant beam radius for diffusely 
reflected or scattered beams. Since we compare the convergence of 
both algorithms for nonspecular regions, this improvement of PPB 
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(a) ours, 16 grid, 20M beams 


(b) ours, 643 grid, 20M beams 


Fig.21. Influence of the size of the covariance grid over the gain in conver- 
gence. A very coarse grid (16° in this example) will conservatively spread 
high-frequency values into regions where the radiance is actually low fre- 
quency, failing to improve the convergence in these regions. Note: we chose 
to use an overly coarse grid in this example to make the effect more apparent. 


was not necessary. Note that both algorithms would benefit equally 
from adding ray differentials. 

We use a sufficiently large starting radius to improve the con- 
vergence of indirect effects while still keeping convergence of the 
direct part. In all our examples, specular paths from the light are 
converged. 


7.2 Limitations 


The various techniques we present, based on our frequency analysis, 
effectively improve convergence in regions of the volume where un- 
dersampling can be performed without loss of accuracy. If frequency 
is high everywhere—such as in a very highly varying medium or 
in a volume where a large number of small shadow rays are cast— 
our a priori analysis naturally predicts that the sampling needs to 
be uniformly dense. In this case, the computation of covariance 
information would naturally not help improve convergence. 

Using volumetric covariance implies an approximation, since it 
neglects the angular covariance of the incident light. Our method 
captures variations in the volumetric fluence that, for reasonably 
nonspecular phase functions, remains close to the variance of the 
radiance while only requiring a small storage cost. In the case of 
a very specular phase function at the last bounce before the camera, 
a 3D covariance grid is likely to produce a conservative overestima- 
tion of the bandwidth. A more accurate approach would require also 
storing directional information in the covariance grid and would not 
invalidate our frequency analysis. 

The size of the covariance grid is another important limitation 
as it determines the scale at which we can optimize the radius 
reduction. A coarse grid will conservatively estimate small kernel 
sizes in low varying regions since high frequencies will leak outside 
of the region where they actually take place (this is illustrated in 
Figure 21). 

The paraxial approximation used by the frequency analysis lim- 
its the capacity of our predictions to describe the full directional 
variations of light with a few photons. The paraxial approximation 
is valid for angles below one degree. However, using our estimates, 
based on this assumption, to estimate light variations works in 
practice. 
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7.3 Discussion 


Performing our analysis in the Fourier domain around light paths 
brings us a local characterization of the signal’s bandwidth. 
Wavelets also bring bandwidth estimation of a signal. However, they 
perform a local analysis at the cost of making operations like con- 
volution and scaling much more complex. In our case, localization 
is already brought by windowing our analysis around a particular 
point of the signal, and the Fourier transform appears to be the most 
simple approach to characterize bandwidth. Polynomial bases, in 
turn, don’t offer simple expressions for convolution and scale. 

Our theoretical derivations perform first-order approximations of 
the scattering and absorption operators, as opposed to first-order 
approximations of the spectra. Linearly approximating the spec- 
tra would be meaningless; in our framework, spectra contain all 
possible frequencies. The only assumption made is the paraxial 
approximation [Durand et al. 2005]. 

As a comparison to the work of Yue et al. [2010], our adap- 
tive sampling strategy is based on the variance of the illumination, 
whereas their algorithm is based on the maximum variance of the 
density in the medium along a light path. Therefore we avoid wast- 
ing time oversampling highly varying regions with low energy. 
While adaptive sampling techniques are usually based on an a pos- 
teriori estimation of the energy (sometimes the bandwidth) of the 
signal, we base our sampling on an a priori prediction of the variance 
of the signal. 

Kulla and Fajardo [2012] propose strategies for importance sam- 
pling participating media. Our approach is complementary since 
we believe that importance sampling metrics can be derived from 
the volumetric covariance. Combining our covariance prediction 
tool with their importance sampling metrics would allow to drive 
the importance sampling by the actual distribution of energy in the 
solution. 

Keeping a large radius for progressive photon beams could slow 
down the selection process if using an acceleration structure (such as 
a KD-tree [Sun et al. 2010]), since this increases the branching factor 
of the KD-tree search. In our CUDA implementation, that follows 
the method described by Jarosz et al., not having an acceleration 
structure removes this penalty. 

When filling the covariance grid, we do not need to record a direc- 
tional distribution of incident illumination to weight the covariance 
contributions from incident rays, since these rays are path-traced 
and arrive with a probability that is already proportional to the 
illumination. Consequently, even in its current and simple imple- 
mentation, the covariance grid allows us to perform anisotropic 
filtering of light beams. This is visible in the caustic scene 
(Figure 8) where we can see that the covariance estimates capture 
the directional discontinuities of the light distribution. 

We chose to apply our framework to accelerate progressive pho- 
ton beams rather than the more recent progressive virtual beam 
lights [Novak et al. 2012a]. The two methods, however, use a sim- 
ilar iterative radius reduction scheme. Therefore, one can expect 
to improve on the latter the same way we improve on progressive 
photon beams, using our radius reduction stopping criterion. 

Similar to Novak et al. [2012a], we could use our method to only 
reconstruct the multiple-scattering component of the illumination 
and not store the specular contribution in the covariance grid. Al- 
though we did not do that, we still produce convincing results as 
our method adapts to any effect (be it direct or indirect). Treating 
indirect scattering independently would enhance the speedup factor 
of our method as the reconstructed signal would be of even lower 
frequency. But this would increase the required engineering for the 
implementation (multiple photon maps would be required). 
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Our improved progressive photon beams method removes the 
need to finely tune the parameters of the radius reduction. This, in 
a way, is similar to the Adaptive Progressive Photon Mapping of 
Kaplanyan and Daschbacher [2013]. One notable difference is that 
our method does not need to evaluate the Hessian of the radiance 
using density estimation and, as shown in Section 6.2, our estimate 
is much more robust. 

Adding the time analysis [Belcour et al. 2013] is straightforward 
but currently limited to rigid motion. The analysis of time varying 
media is also possible, but beyond the scope of this article. Time 
analysis could be implemented using the 3D covariance grid by 
integrating the time dimension. This way, motion events are blurred 
according to the resulting appearance. A 4D grid would be necessary 
to perform temporal coherent filtering. Adding depth of field is also 
orthogonal to this work, but we expect it would not cause particular 
issues. 


8. CONCLUSION 


We proposed the first extension to participating media of the Fourier 
analysis of local light fields. This is a very important problem 
amenable to performance acceleration, since participating media 
typically have lower frequencies. 

We show how to extend the use of covariance matrices in a 
principled manner to represent the spectrum of the light field 
including scattering and absorption. We derive the equations to 
combine the information carried by each light path into a set of 
3D frequency prediction metrics and to compute them from a 
common quantity, namely the volumetric covariance, stored in a 
grid. We used these metrics to improve the convergence and effi- 
ciency of image-space adaptive sampling and reconstruction, cam- 
era ray integration, and for accelerating the progressive photon beam 
approach. 

Several future avenues of research exist. While we demonstrated 
the use of this analysis for the gather part of photon mapping, our fre- 
quency estimate could also be used to drive photon tracing, like, for 
example, using a function of frequency as the acceptance function 
of Metropolis photon sampling [Fan et al. 2005]. Another interest- 
ing research avenue would be to extend our analysis to oriented 
media [Jakob et al. 2010]. 

We did not use an adaptive covariance grid. To do this, the local 
resolution of the covariance grid needs to be adapted to the variations 
of the covariance information to be stored, and would spare the 
need to specify an initial resolution. We leave this question to future 
work. 
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APPENDIX 


A. LAPLACIAN FROM COVARIANCE 
OF SPECTRUM 


Let f be a function defined over a 4D domain. Following basic 
Fourier theory, the DC of the second derivative is the integral of its 
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spectrum 


2 2r 
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where the integral is carried over the entire 4D Fourier domain. We 
expand the Fourier transform of the derivative over each variable: 


af 
Vi, j = (0 
tad bday ) 
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—4r? f Q;2; f(QaAdQ. 


This formulation almost fits the covariance of the power spectrum. 
It actually does when F(Q) is real and positive, and we miss the 
normalization factor (which is f F = f(0) = 1). There exist a wide 
class of such functions [Giraud and Peschanski 2006], the simplest 
of which are Gaussians. In this case we have 


rf 


Ox; 0X; 5 


—4n? f(O) Eixo) = 


In practice, this means that we can approximate the covariance of the 
power spectrum of a function in a small window around a point, as 
soon as the function is close to the osculating Gaussian at that point. 
For 4D functions, the covariance matrix of the windowed power 
spectrum of the function around a particular point (corresponding 
tox — f(X +x)) is therefore well approximated by the following 
diagonal matrix in the frame of the principal curvatures: 


In the most general case, the Hessian matrix must be diagonalised 
before taking the absolute values on the diagonal. Similarly, we 
have 


ə 2 
Dx) © | 4 


Jx? (Xo) 


1 
4n?| f(0)| | 


x Ar? f (Xo)Tr(Z(Xo)). (21) 


B. LOCAL FOURIER COVARIANCE 
OF H-G FUNCTION 


The Henyey-Greenstein phase function is defined by 


1 1- g? 
4T (1 + g? — 2gc)? 


Pg(@j, Ws) = with c = @;.@; = cosa. 


We are interested in the 2D covariance of the windowed spectrum 
of function p, defined by 


P80, db) = pg(cos(a@ + 5A) cos(s@)). 


We recall that the 2D covariance matrix of a phase function p, 
noted o,, is defined by Eq. (13): 


1 
= 4725(0, 0) 


#20, o) 0 
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Taking the second derivatives of p and evaluating them at 60 = 
5 = 0 gives 


ap (0,0) = 3g(g? — 1)(2(g + 1) cosa + g(3 cos(2@) — 7)) 
0602-7 82(g2 — 2g cosa + 1)7/2 
ap 
950050 Č = 8 
ap 3g(g? — 1) cosa 
770,0) = 2 5/2" 
ad 47 (g? — 2g cosa + 1)5/ 


The Hessian being diagonal, we don’t need additional rotations 
before taking the absolute value. When w = 0, the second deriva- 
tives match (to 3g(¢+1)/(42(g—1)*)), as expected for a rotationally 
symmetric function. To get the covariance in Eq. (14), we further 
divide by 417p(0, 0) = 47? p; (a). 
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